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In networks with small-world topology, which are characterized by a high clustering 
coefficient and a short characteristic path length, information can be transmitted efficiently 
and at relatively low costs. The brain is composed of small-world networks, and evolution 
may have optimized brain connectivity for efficient information processing. Despite many 
studies on the impact of topology on information processing in neuronal networks, 
little is known about the development of network topology and the emergence of 
efficient small-world networks. We investigated how a simple growth process that favors 
short-range connections over long-range connections in combination with a synapse 
formation rule that generates homeostasis in post-synaptic firing rates shapes neuronal 
network topology. Interestingly, we found that small-world networks benefited from 
homeostasis by an increase in efficiency, defined as the averaged inverse of the shortest 
paths through the network. Efficiency particularly increased as small-world networks 
approached the desired level of electrical activity. Ultimately, homeostatic small-world 
networks became almost as efficient as random networks. The increase in efficiency 
was caused by the emergent property of the homeostatic growth process that neurons 
started forming more long-range connections, albeit at a low rate, when their electrical 
activity was close to the homeostatic set-point. Although global network topology 
continued to change when neuronal activities were around the homeostatic equilibrium, 
the small-world property of the network was maintained over the entire course of 
development. Our results may help understand how complex systems such as the brain 
could set up an efficient network topology in a self-organizing manner. Insights from our 
work may also lead to novel techniques for constructing large-scale neuronal networks by 
self-organization. 
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1. INTRODUCTION 

The synaptic wiring of cortical networks is key to the function- 
ality of the brain and a precondition for all cognitive behavior 
(Park and Friston, 2013). How are synaptic connections set up 
between the brain's billions of neurons so that cost efficiency 
(Latora and Marchiori, 2001), for example in terms of wiring 
length or energy consumption associated with information trans- 
mission, is maximized? Due to the presence of long-range con- 
nectivity, the brain can be regarded as a highly efficient graph, 
in the sense that information has to pass only a few intermedi- 
ate neurons ("nodes") to travel across the whole brain (Kaiser 
and Hilgetag, 2006). At the same time, brain connectivity is 
local and clustered (Hilgetag and Kaiser, 2004), making net- 
works less vulnerable because of the many alternative routes 
that exist between two nodes. Network topology that combines 
these two properties, long-range connectivity and high cluster- 
ing, is called small-world (Watts and Strogatz, 1998). Small-world 
networks have the advantages of local connectivity combined 
with a high efficiency brought about by a small number of 



long-range connections (Watts and Strogatz, 1998). The brain 
seems to be optimized for maximizing cost efficiency of par- 
allel information processing (Achard and Bullmore, 2007) by 
widely adopting small-world topology (Sporns and Zwi, 2004; 
Bassett and Bullmore, 2006). Already in the infant brain, con- 
nectivity has small-world characteristics (Fransson et al., 2011). 
A recent fMRI study reported an increase in small-worldness of 
brain networks in the first 2 years of life that goes along with a 
growing number of long-distance connections and therefore an 
increase in global efficiency (Gao et al., 2011). So far, analysis of 
topology has focused on neuronal networks with static connec- 
tivity, in which plasticity arises from changes in the strength of 
existing synapses rather than from the rewiring of connectivity. 
However, particularly during development but also in adult- 
hood, connectivity is not fixed (Butz et al., 2009b), and synapse 
formation goes along with massive synapse deletion and reorga- 
nization of connectivity (Missler et al, 1993; Siegel and Lohmann, 
2013). This observation triggered the question of this study: how 
does network topology develop during ontogeny when synapse 
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formation and pruning cause a constant rewiring of network 
connectivity? 

Various explanations have been proposed to account for the 
development of synaptic connectivity. For example, axon guid- 
ance molecules may form the basis of a genetically-encoded 
developmental scheme (Yamamoto et al, 2002; Borisyuk et al., 
2008). Target neurons may secret signaling molecules that can 
attract or repel axons. Axons can then follow or move away from 
the concentration gradient. This form of chemotaxis is usually 
discussed in the context of the formation of global connectiv- 
ity. However, for the formation of local connectivity chemical 
cues are less suitable, since they fail to establish stable gradients 
over very short distances, below 0.7 mm (Kaiser et al, 2009). 
Synaptic adhesion molecules (e.g., neuroligins) were proposed as 
molecular cues for local synapse formation (Scheiffele et al., 2000; 
Stan et al, 2010). In addition, mechanical forces in the tissue 
could influence neurite outgrowth (Franze, 2013). Alternatively, 
the formation of local connectivity may basically be random 
(Braitenberg and Schtiz, 1998), just depending on the acciden- 
tal overlap of axons and dendrites (Binzegger et al., 2004; van 
Pelt and van Ooyen, 2013; McAssey et al., 2014; van Ooyen et al., 
2014). With random synapse formation, the chance of forming 
connections decreases with increasing distance between neurons 
(Kaiser et al, 2009). 

Although random overlap of axons and dendrites may explain 
emerging connectivity, it does not account for the actual driv- 
ing forces underlying neurite outgrowth. There is ample evidence 
that electrical activity shapes neuronal morphology (Dalva et al., 
1994; Wong and Ghosh, 2002; Uesaka et al, 2005; Butz et al., 
2009b) and network formation (Ko et al, 2013). Electrical activ- 
ity influences neurite outgrowth and retraction (McKinney et al., 
1999a; Konur and Ghosh, 2005; Lohmann and Wong, 2005; Tailby 
et al., 2005; Hutchins and Kalil, 2008), as well as the formation 
and deletion of axonal boutons and dendritic spines (McKinney 
et al, 1999b; Groc et al, 2002; Jourdain et al, 2003; Kirov et al, 
2004; Hofer et al., 2009). Experimental findings further suggest 
that neuronal morphogenesis is driven by the need of neurons to 
establish and maintain a homeostatic equilibrium of their aver- 
age electrical activity (Kirov et al, 2004; Keck et al, 2008, 2013). 
Restoration of neuronal firing rate after a change in neuronal 
input has been found experimentally after, for example, focal reti- 
nal lesions (Hengen et al., 2013). Based on these experimental 
findings, we postulated that whenever during development (van 
Ooyen and van Pelt, 1994; Van Ooyen et al, 1995; Tetzlaff et al, 
2010) or in the mature brain (Butz et al, 2008; Butz and van 
Ooyen, 2013) a neuron senses a deviation of its electrical activ- 
ity from a homeostatic set-point, it will initiate changes in its 
morphology that increase the chance of synapse formation or 
break existing connections so that its firing rate may be restored. 
Here we investigate what the impact is of neurons regulating their 
electrical activities homeostatically on network formation and 
emerging network topologies. 

In order to study the impact of electrical activity on emerg- 
ing network topology, we used our recent Model of Structural 
Plasticity (MSP) (Butz and van Ooyen, 2013). There are impor- 
tant earlier models of homeostatic structural plasticity, such as 
the compensation model by Dammasch et al. (Dammasch et al., 



1986, 1988; Cromme and Dammasch, 1989; Butz and Teuchert- 
Noodt, 2006; Butz et al, 2006) and the activity-dependent neurite 
outgrowth model by van Ooyen (van Ooyen and van Pelt, 1994; 
Van Ooyen et al., 1995). The latter model, which studied the 
reciprocal interactions between neuronal activity and network 
formation, successfully accounted for experimental data on devel- 
oping cell cultures (van Ooyen and van Pelt, 1994; van Oss and 
van Ooyen, 1997; Abbott and Rohrkemper, 2007; Tetzlaff et al, 
2010). However, both models are not suited for studying topol- 
ogy development. The compensation model lacks topology at 
all, whereas the representation of neuritic fields by circles in the 
model by van Ooyen imposes too strong constraints on net- 
work topology. In van Ooyen's model, neurons always connect to 
their direct neighbors before connecting to more distant neurons. 
Therefore, we developed our MSP (Butz and van Ooyen, 2013), 
in which we replaced the circle representation by discrete synaptic 
elements whose numbers change in an activity-dependent man- 
ner. Synapses are formed in a random and distance-dependent 
way by combining synaptic elements from different neurons. 

In MSP, the local activity-dependent growth process in com- 
bination with a simple kernel function favoring the formation 
of short-range connections over long-range connections shapes 
the development of small-world networks. Our simulation results 
revealed an interesting property of homeostatic growth: as soon 
as most neurons approached homeostasis in electrical activity, 
they started forming more long-range connections than expected 
from the kernel function. Although the clustering coefficient 
decreased as a result of the formation of long-range connections, 
the network maintained its small-world property. Furthermore, 
connectivity became more diverse, as indicated by a decreasing 
betweenness centrality, and attained a higher global efficiency 
(defined as the averaged inverse of the shortest paths between 
all neurons in the network) than small-world networks with- 
out homeostasis. Our findings may account for experimental 
data on the topology of developing dissociated cell cultures 
(Downes et al., 2012). Interesting similarities were also found 
between the model and early human brain development (Gao 
et al., 2011) with respect to an increasing number of long- 
range connections and an increasing global efficiency during 
development. 

2. MATERIALS AND METHODS 
2.1. NEURON MODEL 

The model network consisted of n = 400 neurons, of which 80% 
were excitatory and 20% inhibitory. Excitatory neurons were 
placed with some jitter on a 20 x 16 grid with a spatial dis- 
tance between two grid points of 150 fim. Inhibitory neurons 
were placed evenly between the excitatory neurons. The same net- 
work layout was used as in Butz and van Ooyen (2013). We used 
Izhikevich's model (Izhikevich, 2003) to simulate neuronal elec- 
trical activity. This model has two differential equations, one for 
the membrane potential v (in mV) and one for a recovery variable 
u (in mVms -1 ), enabling re-polarization after an action potential: 

* = k lV 2 + k 2 v + k 3 -u + I 

%=a(bv-u) (1) 
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where k\ = 0.04mV _1 ms , k 2 = 5 ms , /c 3 = 140mVms~ 
and t is in ms. Every time a neuron fires (v > 30 mV), v and u 
are reset: 



if v > 30 mV, then 



f v <- c 
( u <— u 



+ 



(2) 



As in Izhikevich (2003), the following parameter values were used: 
a = 0.1ms" 1 , b = 0.2 ms" 1 , c = -65 mV, d = 2mVms~ 1 . We 
used the same dynamics for excitatory and inhibitory neurons. 
Each neuron receives input I = I syn + I ext , consisting of synap- 
tic input Ism from within the network and external input I ext . 
Neurons interchange electrical signals on a millisecond timescale 
without synaptic delay. Synaptic input consists of the incoming 
action potentials from the presynaptic neurons low-pass filtered 

by an exponential filter function h(t) = exp ^— j^j with decay 
constant fj, = 5 ms. Network connectivity Wu is defined as the 
number of synapses from neuron j to i. If a synapse exists, it 
has a fixed strength of 1 mVms - 1 . Neurons are either excitatory 
or inhibitory. Indices refer to excitatory neurons if i or j e {Ex} 
and to inhibitory neurons if i or j e {In}. As in Izhikevich (2003) 
and Butz and van Ooyen (2013), I ext is delivered as white noise 



with mean 5 mVms and standard deviation 1 mVms 



In 



some cases lower input values were required and specified where 
relevant. 

The intracellular calcium concentration of a neuron is used as 
a low-passed filtered average of its firing frequency (Butz and van 
Ooyen, 2013). Every time a neuron fires, the calcium concentra- 
tion is increased by a fixed amount; otherwise the concentration 
decreases exponentially to zero: 



d[Ca 2 +]. 
dt 



[c* 2+ ], 



+ P if v > 30mV 



(3) 



otherwise 



where /3 = 0.001 ms" 1 and rc„ = 10,000 ms. We defined a set- 
point 6 = 0.7 in calcium concentration, corresponding to an 
intermediate level of average electrical activity. Every time the 
neuron's calcium concentration deviates from e, it will induce 
structural changes in connectivity to restore the desired average 
level of activity, as described below. 

2.2. SYNAPSE MODEL FOR STRUCTURAL PLASTICITY 

During development, neurons show a pronounced formation 
and pruning of synapses. To simulate reorganization of synaptic 
connectivity, standard models of synaptic plasticity, in which con- 
nectivity is considered fixed with plasticity merely arising from 
changes in the strength of existing synapses (modeled as weight 
factors), are not suitable. For this study on the development of 
neuronal networks, we therefore used our Model of Structural 
Plasticity (MSP) (Butz and van Ooyen, 2013). The characteristic 
feature of this model is that it represents synapses as consisting of 
two recombinable synaptic elements, namely an axonal element 
and a dendritic element. Axonal elements represent countable 
presynaptic specializations for transmitter release such as axonal 
terminals or boutons, while dendritic elements represent the 



postsynaptic counterparts, i.e., postsynaptic receptor plates on 
dendritic spines or the dendritic shaft. Axonal and dendritic ele- 
ments are either excitatory or inhibitory. Any neuron i can form 
A,- axonal elements, which are excitatory if the neuron is excita- 
tory, or inhibitory if the neuron is inhibitory. At the same time, 
a neuron, irrespective of its type, can express Df excitatory and 
Df inhibitory dendritic elements. Complementary elements can 
merge to form a synapse (excitatory axonal with excitatory den- 
dritic elements, and inhibitory axonal with inhibitory dendritic 
elements). Synaptic elements form and delete independently from 
a synaptic contact partner. In case a neuron deletes a synaptic 
element that is bound in a synapse, the complementary synap- 
tic element on the other neuron remains and becomes vacant and 
available again for synapse formation with a new target. Therefore 
MSP allows for synaptic rewiring. 

2.3. HOMEOSTATIC GROWTH RULES 

It is well documented that neurons change their morphology in 
an activity-dependent fashion during development (Butz et al., 
2009b; van Ooyen, 201 1). The way in which neurons change their 
morphology suggests that they try to maintain homeostasis of 
electrical activity (van Ooyen and van Pelt, 1994; Van Ooyen 
et al, 1995; Butz et al, 2009a; van Ooyen, 2011; Butz and van 
Ooyen, 2013). That is, neurons in which the average activity is 
too low start forming new neuritic structures, whereas neurite 
growth is halted or neurites are pruned when activity is higher 
than a desired level (homeostatic set-point). In our MSP, we 
abstracted away from describing detailed neuronal morphology 
and assumed that changes in morphology effectively result in 
changing numbers of axonal and dendritic elements — the con- 
tact sites for synaptic connections. Homeostatic adaptation of the 
number of axonal and dendritic elements was modeled by the fol- 
lowing growth rule (Figure 1), in which dz represents the change 
in the number of synaptic elements, with z being A, D ex or D 1 ": 



dz 
dt 



([c fl 2 +H)/<u 



l + e 



(4) 



where e is the set-point of the intracellular calcium concentra- 
tion, corresponding to a desired average firing rate of the neuron, 
and v is the growth rate of synaptic elements. We chose v = 
10 _4 ms _1 , which is slow enough that electrical dynamics and 
structural dynamics do not interfere, yet fast enough not to slow 
down the simulations unnecessarily. For reasons of simplicity and 
because of a lack of detailed experimental data, we applied iden- 
tical sigmoid growth rules to all types of synaptic elements. In 
a random and distance-dependent recombination process, newly 
formed synaptic elements were distributed to matching synaptic 
elements, thereby forming synapses and creating the pattern of 
connectivity in the network. 

2.4. KERNEL FUNCTION FOR SYNAPSE FORMATION 

As in our previous work on MSP (Butz and van Ooyen, 2013), 
we assumed that synapse formation is more likely between adja- 
cent neurons than between distant ones. We applied a two- 
dimensional Gaussian-shaped kernel function centered at the 
x, y-coordinates of neuron z, with Ku the distance-dependent 



Frontiers in Synaptic Neuroscience 



www.frontiersin.org 



April 2014 | Volume 6 | Article 7 | 3 



Butz etal. 



Homeostatic small-world networks 



1 














dz 
df 








0 














-1 


■ 












0 0.5 £ 1.0 [Ca 2+ ] 1-5 


FIGURE 1 | Identical sigmoidal growth rules were used for all types of 
synaptic elements to determine the change dz/dt in the number of 
synaptic elements in dependence on the intracellular calcium 
concentration ^Ca 2+ J. z needs to be replaced by the respective type of 
element A, D ex or D ln . The homeostatic set-point e is the value of the 
calcium concentration where dz/dt = 0. 



likelihood for synapse formation between neuron j and 



(pos xj - pos xi ) 2 + (posyj - posyj) 2 



K, 



(5) 



where pos x i is the x-coordinate and pos yi is the y-coordinate of 
postsynaptic neuron i, and pos x j and pos y j are the coordinates 
of presynaptic neuron j. The probability for autapse connections 
(i.e., a neuron connecting to itself) was set to zero (Ku = 0 for 
i= j). For these simulations we chose a = 1 X 150 /xm where 
150 /im is the distance between two grid points. 

In order to investigate the impact of this distance-dependency 
on emerging network topologies, we additionally grew networks 



with a flat kernel, i.e., with K, 



'J. >¥^j 



1 and K, 



0 for i 



2.5. SYNAPSE FORMATION AND DELETION 

The MSP algorithm (Butz and van Ooyen, 2013) proceeded 
in three steps to update network connectivity in an activity- 
dependent fashion. First, electrical activity of all neurons was 
computed continuously over time. Secondly, depending on the 
average electrical activity and according to the homeostatic 
growth rule (Equation 4), the number of elements changed con- 
tinuously, too. Thirdly, at discrete time points, network connec- 
tivity W was updated by synapse formation and deletion. Because 
of the low growth rate v = 10~ 4 ms -1 , changes in connectivity 
are very slow compared to changes in activity, so it was not nec- 
essary to update connectivity at every time step but only at every 
100 ms. The timescale of network formation in the model corre- 
sponds in principle to a timescale of days or weeks. However, as 
described above under homeostatic growth rules, the timescale of 
structural changes was chosen so that the simulations were not 
unnecessarily slowed down, which means that the total duration 



of the simulation in milliseconds is not required to sum up to 
weeks. 

2.5.1. Synapse deletion 

Since network connectivity is updated at discrete time steps but 
synaptic elements change continuously over time due to the 
activity-dependent growth rules, it can happen that a neuron has 
more outgoing synapses than axonal elements or more incoming 
synapses than dendritic elements at the time of the next update 
in network connectivity. In this case, the neuron has to delete the 
surplus of synapses and to update connectivity. 

To update connectivity, the algorithm needs to select which 
synapses are to be removed. All synapses have an equal chance of 
being deleted. Note, however, that multiple synapses can co-exist 
from neuron j to i and that the more synapses there are, the higher 
the chance that a synapse between neuron j and i will be deleted. 
The probability Pfj for synapse deletion between neuron and i 
is computed by the following master equation that captures four 
different cases: 



ndel 



(6) 



For deletion of incoming synapses, we need to distinguish 
between excitatory and inhibitory synapses in Equation (6). For 
deleting incoming excitatory synapses of neuron i e {In U Ex}, we 
sum up Wk,j over all I e {Ex}. For deleting incoming inhibitory 
synapses of neuron i e {In U Ex}, we sum up W^j over all I e {In}. 
For deletion of outgoing excitatory synapses of excitatory presy- 
naptic neuron with ; e {Ex}, in Equation (6) all synapses are 
considered to any postsynaptic neuron fc with fe € {In U Ex}. 
Thus, we sum up W^i over all k e {In U Ex}. The same holds true 
for outgoing inhibitory synapses with j e {In}. 

Sequentially, outgoing and incoming excitatory and inhibitory 
synapses were selected for deletion. For every type of synapse, the 
accumulated sum of Pfj [see description of Equation (6) for the 
range of i and j] gave a probability distribution from which we 
drew the required number of synapses to be deleted. The selected 
synapse was deleted by reducing the respective entry W;j in the 
connectivity matrix by one. It can happen that more then one 
synapse is selected for deletion from the same connection; to i. In 
this case, the implementation of the algorithm makes sure that the 
number of synapses to be deleted did not exceed W,- 1. Whenever a 
neuron deletes a synaptic element that is bound in a synapse, the 
complementary synaptic element on the other neuron remains 
and becomes vacant again. 

2.5.2. Synapse formation 

For synapse formation, the algorithm checked whether a neuron 
gained vacant synaptic elements, i.e., whether the total number of 
synaptic elements exceeded the number of bound synaptic ele- 
ments of this type. Matching vacant synaptic elements (vacant 
excitatory axonal elements A^' vac with vacant excitatory dendritic 
elements T).' 9C , and inhibitory axonal with inhibitory dendritic 
elements rj m - vflc ) were randomly distributed among each other 
with probability density function pf orm . The probability F°™ 
for forming new synapses between neuron j and i depended on 
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the number of vacant synaptic elements they offered and on the 
Euclidean distance between neuron j and i: 



pform 



j€{Ex) : 



/€ {IB} : 



j^vac j^ex.vac 



A vac \ 1 p,ex,vac -^y 



A vac r-iin,vac 9 



withie {Bell In}. 



(7) 



The minor number of vacant excitatory and inhibitory axonal 
or dendritic elements determined how many new excitatory 
and inhibitory synapses, respectively, could at most be formed 
(so-called potential synapses) in every update of connectivity. 
Thus, the number of excitatory and inhibitory potential synapses 
equaled 



^j-PotSyn,ex _ 



^j-PotSyn, in 



1 (l2ie{Ex 



€{ExUIn) k 



i I A vac rtin,vac \ 

' yZ^ie{In}< A i ' l^K(E{ExUIn) u k ) 



(8) 



for every update in connectivity. 

From this distribution, the algorithm chose at maximum 

M PotSyn,ex excitatory an d M PotSyn,in inhibitory connections at 

which new synapses were created. The respective entries Wij in 
the connectivity matrix were then increased by one. A connection 
was chosen by drawing a random number from a uniform dis- 
tribution and comparing it to the accumulated probabilities p[° rm 
for all excitatory connections and all inhibitory connections of the 
entire network. That connection was chosen that had the high- 
est accumulated probability that the random number just did not 
exceed. If, for this try, the random number exceeded all accumu- 
lated probabilities, no synapse was formed. Hence, not necessarily 
all of the potential synapses were formed. 

Additionally, synapse formation needed to fulfil the condition 
that the number of newly-formed synapses from neuron j 
to i did not exceed the number of vacant synaptic elements that 
neuron; and i offered: 



j e {Ex} 



j e {In} 



mm(AJ ac 



min(A vflc 



jyn,vac\ 



with % 6 {ExU In}. 



(9) 

In every update, this condition was checked and synapse for- 
mation infringing this condition was rejected. Alternatively, the 
update of connectivity can also be implemented in a purely local 
fashion (Butz and van Ooyen, 2013). For small networks, the 
current implementation is more efficient than the original MSP 
algorithm (Butz and van Ooyen, 2013) because run time is not 
dependent on the growing numbers of vacant synaptic elements 
but proportional to the square of the number of neurons in 
the network. However, for large networks with n >> 1000, n 2 
will quickly out-range the number of vacant synaptic elements, 
in which case looping over the number of synaptic elements 
would be faster. Particularly in Matlab, the current description of 
MSP allows for an elegant vectorized implementation providing 
additional speed up. 



2.6. DEVELOPMENT OF NON-HOMEOSTATIC NETWORKS 

To investigate the impact of homeostasis in electrical activity 
on developing network topology, we created for every home- 
ostatic network a corresponding non-homeostatic network. At 
every update in connectivity, we took the number of synapses 
from the homeostatic network and distributed them in the non- 
homeostatic network with the same kernel function as used in the 
homeostatic network. Hence, the placement of synapses in the 
non-homeostatic network was purely dependent on the kernel 
function but did not meet the homeostasis criterion. The algo- 
rithmic implementation for placing synapses was the same for 
homeostatic and non-homeostatic networks, with PvJ" 1 = Kj j. 

Instead of distributing M FotSyn synapses (Equation 8), we simply 
distributed the total number of synapses from the homeostatic 
network. Since synapse formation was not limited by numbers of 
vacant elements, Equation (9) was not applicable. 

2.7. TOPOLOGY MEASUREMENTS 

A neuronal network can be seen as a graph with the neu- 
rons being the nodes and the synapses being the edges or links 
between two nodes. Since the presynaptic neuron always activates 
the postsynaptic neuron (and never the other way around), we 
regard the graph as directed. At every update in connectivity, we 
assessed those graph theoretical measures that are indicative of 
small-worldness and network efficiency. In addition, betweenness 
centrality was measured to determine the importance of nodes 
in the network. To reduce the complexity of the assessment, we 
considered only the topology of excitatory synaptic connections 
W ex ' ex between the n ex excitatory neurons. For the graph theo- 
retical assessments, the brain connectivity toolbox by Olaf Sporns 
et al. (Rubinov and Sporns, 2010) was used. 

2.7. 1. Weighted characteristic path length 

The characteristic path length L is the average shortest path from 
one node to any other node in the network. Path length is defined 
as the number of links that need to be traveled to go from one 
node (possibly via intermediate nodes) to any other node. On top 
of this definition, a direct link between two nodes in a weighted 
network is considered "shorter" the stronger the weight of the 
link is. For our network, we take the number of synapses Wj X ' ex 
between two directly linked neurons; and i, with i, j e {Ex}, as the 
weight of the connection and the inverse \/W^' ex as the length Z,- j 
of the connection. The shortest path d u] is then the smallest sum 
of connection lengths that lead from neuron; to i via any interme- 
diate neurons. Our calculation of the weighted characteristic path 
length was based on an implementation of Dijkstra's algorithm 
for computing the shortest path in weighted directed networks by 
Rubinov and Sporns (2010). 

2.7.2. Weighted clustering coefficient 

The clustering coefficient is an indication of how strongly nodes 
in a network are interconnected. It can be measured by the num- 
ber of triangles, t D , one node forms with any other two nodes 
in the network divided by the total number of possible triangles, 
TP. For weighted and directed graphs, one needs to realize that 
the adjacency matrix (in our case, the connectivity matrix W ex ' ex 
of the excitatory neurons) is not symmetric and that the entries 
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of the adjacency matrix are not one but can have any weight. 
According to Fagiolo (2007), the clustering coefficient of a single 
node in a weighted directed network Cf is computed as 



Cf (W ex ' ex ) 



_i 

jD 



+ 



f^ex.ex^A 



1/3]' 



2[df 



I) -2d?} 



(10) 

the k th root of the entries of the 



where (W ra - ra )[V« 

matrix for i,je {Ex}, and (W ex,ex ) T is the transposed W ex ' ex 
matrix. The variable d m denotes the total degree of node i 
(the degree counts the number of either incoming or outgo- 
ing edges per node, and the total degree is the sum over both 
the incoming and outgoing edges), and d** stands for the num- 
ber of bilateral edges of node i (the number of nodes node i 
projects to and receives edges from, excluding self-interactions of 
node i). 

The overall clustering coefficient of the network is thus C D = 
j C?. Note that for this assessment we only con- 
sidered n ex excitatory nodes and excitatory connections W ex ' ex . 
For a more detailed description of Equation (10), see Fagiolo 
(2007). We computed the clustering coefficient of the develop- 
ing neuronal networks at every update in connectivity by using 
the implementation by Mika Rubinov from the brain connectivity 
toolbox (Rubinov and Sporns, 2010). 

2.7.3. Small-world parameter 

To estimate small-worldness of the developing networks, we 
applied the formalism by Humphries and Gurney (2008): 
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We replaced the clustering coefficient C and the characteristic 
path length L by the corresponding versions for weighted directed 
graphs as described above. C mnd and L mnd were taken from an 
Erdos-Renyi random graph generated with the same number 
of synapses as in the developing networks at every update in 
connectivity. 

2.7.4. Betweenness centrality 

The local betweenness centrality is a measure for the impor- 
tance of a node in a network. A high betweenness central- 
ity of a node i means that many shortest paths between 
any two nodes k and / pass through node i. Thus, the local 
betweenness centrality counts the number of times, ffki(i), that 
node i is on a shortest path between two nodes k and I. 
The local betweenness centrality is normalized by the num- 
ber (Tj-; of alternative shortest paths between k and / that do 
not pass through node i. Global betweenness centrality is the 
sum of all local betweenness centrality values of the individual 
nodes: 
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(12) 



Consequently, betweenness centrality provides a measure for 
how well networks are interconnected. A high local or global 
betweenness centrality means that individual nodes or the entire 
network, respectively, is badly interconnected, because all infor- 
mation has to travel through the same nodes in the absence 
of alternative routes or by-passes. Clinical data shows that 
after brain lesions, betweenness centrality of directly and indi- 
rectly affected brain areas changes (Wang et al, 2010). Note 
that as for the measurements above, the shortest paths are 
based on weighted excitatory connections W^ x ' ex . Therefore, 
global betweenness centrality was computed by the formalism 
for weighted directed networks by Brandes (2001) as imple- 
mented in the brain connectivity toolbox (Rubinov and Sporns, 
2010). 

2.7.5. Efficiency 

Global efficiency is related to the inverse characteristic path length 
with the advantage over characteristic path length that it can be 
meaningfully computed also of disconnected graphs (Latora and 
Marchiori, 2001; Achard and Bullmore, 2007). While path lengths 
between disconnected cells are infinite, efficiency becomes zero 
and, thus, adds neutrally to global efficiency. 
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(13) 



where n ex is the number of excitatory neurons. 
2.7.6. Euclidean distance 

Although not a topology measure in a strict graph theoretical 
sense, the average Euclidean distance between nodes that are con- 
nected by synapses was measured in order to help interpret the 
development of characteristic path length and clustering coeffi- 
cient. To obtain the average Euclidean distance ED, we multiplied 
the Euclidean distances between all pairs of excitatory neurons in 
the network with the number of synapses between these neurons. 
We then summed up all of the so-weighted distances and divided 
the sum by the number of excitatoy synapses: 



(14) 



with i, j € {Ex} 



3. RESULTS 

We started each network simulation with zero connectivity and 
zero synaptic elements. Due to the homeostatic formation of 
axonal and dendritic elements, model neurons are able to form 
synapses and to adapt the number of synapses so as to reach 
a homeostatic equilibrium in electrical activity (Figure 2). For 
a wide range of set-points e, model neurons adapt their aver- 
age firing rate so that, at the end of each simulation (here 
after 15,000 updates in connectivity), their calcium concentra- 
tion [Ca 2+ ] has converged to e, with a corresponding firing 
rate y (in Hz) that follows the linear relation y = 100 x [Ca 2+ ] 
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(Figure 3). We investigated how the topology of these self- 
organizing networks developed over time when neurons favor 
short-range connections over long-range connections, or, alter- 
natively, when all connections are equally likely. In the first case, 
as we will show in detail later on, networks developed a distinct 
small-world property, with small-world parameter s markedly 
greater than 1, whereas in the second case, s reaches 1. Therefore, 
we will call networks that resulted from distance-dependent 
synapse formation small-world networks and those that resulted 
from synapse formation without distance-dependency random 
networks. 

In networks that favor long-range over short-range connec- 
tions (small-world networks) (Figure 4A), s constantly increased 
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FIGURE 2 | Development of intracellular calcium concentration [Ca 2+ J 
over time. (A) Development of calcium concentration in small-world 
networks arising from synapse formation that favors short-range over 
long-range connections. Mean calcium concentrations averaged over five 
simulations reach the set-point e = 0.7 when the number of synaptic 
elements is regulated homeostatically (yellow), whereas calcium 
concentrations remain much lower when there is no homeostatic 
regulation (gray). (B) Development of calcium concentration in random 
networks without any distance-dependency in synapse formation. The 
random network with homeostatic regulation of synaptic elements (blue) 
also reaches the homeostatic set-point e, whereas the network without 
homeostasis (gray) and the same number of synapses as the homeostatic 
network has much lower values in calcium. Shadow around the curves 
(hardly visible since so small) indicates the standard deviation. 
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Calcium concentration 

FIGURE 3 | Neurons develop their connectivity in order to reach a 
homeostatic set-point e of intracellular calcium. The figure shows the 
firing rates attained for different set-point values of calcium. For each 
calcium concentration, the firing rates of all neurons were pooled from four 
different simulations recorded over the last 20,000 ms. The central mark of 
each box is the median firing rate; boxes represent the 25th and 75th 
percentiles; the whiskers extend to the most extreme data points not 
considered "outliers"; and "outliers" are plotted individually (which show up 
here as thick bold lines). An "outlier" is a value that is more than 1.5 times 
the interquartile range away from the top or bottom of the box. The firing 
rate is proportional to e by a factor of 100. For this set of simulations lower 
external inputs l ext with mean 2 mVms -1 were used. 



and reached a maximum markedly greater than 10 at around 7000 
updates in connectivity. Small-world networks that were set-up 
by an additional rule for homeostasis in electrical activity reached 
a plateau of about s = 10 very early but began to decrease again 
around the time that neuronal activities reached the homeo- 
static set-point e. Nevertheless, small-world networks with home- 
ostasis maintained their small-world property throughout the 
whole course of the simulation (s > 5 at T = 15,000). In net- 
works without distance-dependency in synapse formation (ran- 
dom networks) (Figure 4B), s equaled 1 from the very beginning 
of network development. Random networks with and with- 
out homeostasis did not differ in the course of s (Figure 4B). 
Hence, homeostasis did not seem to have an impact on the 
development of topology in random networks. Knowing that 
homeostasis influenced the development of small-worldness, we 
further analyzed how homeostasis exerted its influence during 
network formation. For some simulations, s could not be com- 
puted for the first few updates in connectivity because of division 
by zero. 

3.1. HOMEOSTASIS INFLUENCES THE CLUSTERING COEFFICIENT AND 
CHARACTERISTIC PATH LENGTH OF DEVELOPING SMALL-WORLD 
NETWORKS 

The mean or characteristic path length in small-world net- 
works without homeostasis showed, after an initial sharp rise, 
a steady decrease and converged toward values of around 3 
(Figure 5A). Homeostatic small-world networks also started 
with a sharp rise and a subsequent decrease in characteris- 
tic path length. The decrease was even steeper than in the 
non-homeostatic case, and the characteristic path length con- 
verged toward slightly lower values than in non-homeostatic 
networks. In addition, networks with homeostasis showed a 
second but minor decrease in characteristic path length when 
activities reached the homeostatic set-point e. The final values 
of the characteristic path length after 15,000 updates in con- 
nectivity were stable in both homeostatic and non-homeostatic 
small-world networks. The initial sharp increase is caused by 
the limited number of synapses in the network at the begin- 
ning of the simulation. Characteristic path lengths in developing 
random networks showed an identical course with and without 
homeostasis (Figure 5B). The final values in random networks 
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FIGURE 4 | Small-worldness of developing networks. (A) Small-world 
networks with (orange) and without (gray) homeostasis in electrical activity. 
(B) Random networks with (blue) and without (gray) homeostasis in 
electrical activity. Means over five simulations per scenario. Shadings of the 
curves indicate standard deviations. 
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were marginally lower than those in homeostatic small-world 
networks. 

The clustering coefficient in developing homeostatic small- 
world networks took a considerably different course from 
the coefficient in small-world networks without homeostasis 
(Figure 6A). Whereas the clustering coefficient in networks with- 
out homeostasis converged, with a small overshoot, toward high 
levels of over 1.6, networks with homeostasis generated a con- 
siderable overshoot and ended up at much lower values as com- 
pared with non-homeostatic networks. After a ramp-up phase, 
homeostatic networks reached a maximum clustering coeffi- 
cient of about one; thereafter, the clustering coefficient decreased 
again with a decreasing negative slope. The maximum cluster- 
ing coefficient was reached when the average electrical activities 
approached the homeostatic set-point e. By contrast, we did not 
see an overshoot in clustering coefficient in homeostatic and 
non-homeostatic random networks (Figure 6B). Therefore, the 
homeostasis in electrical activity had no effect on clustering in 
random networks. 

3.2. HOMEOSTASIS FAVORS LONG-RANGE CONNECTIONS IN 
SMALL-WORLD NETWORKS 

Particularly the development of the clustering coefficient, with its 
pronounced overshoot, determines the emerging small-worldness 
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FIGURE 5 | Characteristic path length in developing networks. (A) 

Small-world networks with (orange) and without (gray) homeostasis in 
electrical activity. The inset in (A) is a close-up of the time interval from 
5000 to 10,000 updates in connectivity that clearly shows the decay in 
characteristic path length in homeostatic networks compared to 
non-homeostatic networks. (B) Random networks with (blue) and without 
(gray) homeostasis in electrical activity. Means over five simulations per 
scenario. Shadings of the curves indicate standard deviations. 
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FIGURE 6 | Clustering coefficient in developing networks. (A) 

Small-world networks with (orange) and without (gray) homeostasis in 
electrical activity. (B) Random networks with (blue) and without (gray) 
homeostasis in electrical activity. Means over five simulations per scenario. 
Shadings of the curves indicate standard deviations. 



in networks favoring short-range connections over long-range 
connections. We therefore further studied how homeostasis influ- 
enced clustering in small-world networks. We tested the hypothe- 
sis that homeostasis produced more long-range connections than 
expected from the kernel function K (Equation 5). Computing 
the average Euclidean distance between the pre- and postsy- 
naptic neuron for every synapse indeed revealed longer average 
Euclidean distances in homeostatic than in non-homeostatic 
networks (Figure 7A). In non-homeostatic networks, the aver- 
age Euclidean distance was constant, because it is directly 
derived from the kernel function (in Equation 5: a = 150 /xm). 
In homeostatic networks, however, we observed two different 
phases. First, during initial network development ([Ca 2+ ] << e), 
the average Euclidean distance converged quickly toward a stable 
value of around 2, which was only slightly higher than in non- 
homeostatic networks. Secondly, when calcium concentrations 
approached the homeostatic set-point e, the average Euclidean 
distances ramped up and reached values greater than 4. The con- 
siderable increase in the average Euclidean distance of synaptic 
connections coincided with a drop in clustering coefficient. As 
the initial high clustering of the network is due to the kernel 
function for synapse formation favoring short- over long-range 
connections, we may conclude that increasing Euclidean lengths 
of synaptic connections give rise to the decrease in clustering 
right at the time neurons approach the homeostatic set-point in 
electrical activity. The effect was also noticeable in the course of 
the characteristic path length but much less pronounced. We did 
not observe a comparable effect of homeostasis on Euclidean dis- 
tances of synaptic connections in random networks (Figure 7B). 

Are the increasing average Euclidean distances of synaptic con- 
nections caused by the fact that average neuronal activities are 
reaching a homeostatic equilibrium, or is this just some net- 
work effect that merely coincides with neurons equilibrating their 
activities? To answer this question, we assessed the number of 
all types of vacant elements (i.e., excitatory or inhibitory axonal 
elements, excitatory or inhibitory dendritic elements) and the 
spatial position of their hosting neurons at every time step of 
the simulation. We first checked whether there was any bias in 
the spatial position of synaptic elements that could generate more 
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FIGURE 7 | Average Euclidean distance between pre- and postsynaptic 
neurons for every synapse during network development. (A) 

Small-world networks with (orange) and without (gray) homeostasis in 
electrical activity. (B) Random networks with (blue) and without (gray) 
homeostasis in electrical activity. Means over five simulations per scenario. 
Shadings of the curves indicate standard deviations. Inset in (A) depicts the 
change in Euclidean distances for homeostatic random (blue) and 
homeostatic small-world networks (orange) with torus boundary conditions. 
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distant connections, e.g., a placement of synaptic elements at the 
boundaries of the network. We accumulated the number of vacant 
axonal (Figure 8A) and dendritic elements (Figure 8B) for each 
neuron for the first 5000 updates in connectivity with [Ca 2+ ] < 
€ as well as for the next 5000 time steps with [Ca 2+ ] ~ e. In 
the beginning of the simulation, we indeed found a little more 
vacant synaptic elements at the network boundaries, which can be 
explained by the fact that neurons at boundaries have less neigh- 
bors than neurons in the center and compensate for this by getting 
a higher number of vacant dendritic elements. However, by the 
time that the homeostatic equilibrium is reached, vacant synaptic 
elements were equally distributed over the network, and therefore 
the placement of vacant synaptic elements cannot account for the 
increasing distances of synaptic connections once the network has 
reached the homeostatic set-point. We found a comparable course 
of Euclidean distances in homeostatic small-world networks with 
torus boundary conditions (see inset of Figure 7A). 

Since we could exclude a spatial bias in the distribution of 
vacant synaptic elements, we tested whether the increasing dis- 
tances of synaptic connections were a direct consequence of the 
activity-dependent growth rules. On the basis of the number of 
vacant elements per neuron and the Euclidean distance between 
the host neurons, we determined the most likely synaptic connec- 
tion every time a new synapses was formed, which is equivalent 
to the maximum of P^™ (Equation 7). Since the change in the 
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FIGURE 8 | Spatial position of vacant synaptic elements in small-world 
networks. The purpose of this figure is to rule out that the position of 
vacant synaptic elements alone gave rise to longer-ranged synaptic 
connections. The left column shows vacant synaptic elements accumulated 
over the first 5000 updates in connectivity, whereas the right column 
shows the accumulation of vacant synaptic elements over the following 
5000 updates. Although some tendency of vacant synaptic elements being 
located at the borders of the network is visible, this effect is gone after 
5000 updates in connectivity. So there is no bias in the distribution of 
vacant synaptic elements when activity reaches the homeostatic set-point, 
and therefore the position of vacant synaptic elements alone cannot 
account for the formation of longer-ranged synaptic connections for 
T > 5000. (A) Vacant axonal elements. (B) Vacant dendritic elements. Color 
scale indicates number of vacant synaptic elements per neuron. 



number of synaptic elements and therefore also the distribu- 
tion of vacant synaptic elements are activity-dependent, the most 
likely synaptic connection to be formed is a direct consequence of 
the neurons' electrical activities. It turned out that in the begin- 
ning of development, when all neurons offer vacant synaptic 
elements, the most likely synaptic connections are those between 
adjacent neurons (Figures 9A,B). In other words, when neuronal 
activities were much lower than the homeostatic set-point e, the 
kernel function had a large impact on the Euclidean length of 
synaptic connections. Therefore, at this stage of development, 
homeostatic and non-homeostatic networks did not differ much 
in connection lengths. However, when the activity of all individ- 
ual neurons approached the set-point 6, vacant synaptic elements 
became rare and matching synaptic elements were available, if at 
all, only between more distant neurons. Expected distances for 
most likely synaptic connections therefore became much larger. 
Moreover, the distribution of expected Euclidean lengths of new 
synaptic connections did not follow the Gaussian-shaped kernel 
any longer but became much wider and flatter. At the same time, 
it took much longer for a synapse to form because, although 
longer-range synaptic connections were the most likely ones, 
their probability was still very low due to the kernel function. 
Due to a lack of shorter-range alternatives, these longer-range 
synaptic connections were nonetheless formed at some point in 
time, because the kernel function is non-zero for all distances. 
Taken together, the Euclidean length of synaptic connections in 
small-world networks is influenced by the homeostatic forma- 
tion of synaptic elements. As expected from the topology data, 
no activity-dependent effect on the Euclidean distance of synaptic 
connections was observed in random networks. 

Additionally, synapse deletion could in principle also influ- 
ence the Euclidean length of synaptic connections if, for example, 
with increasing neuronal activities preferentially short connec- 
tions would be deleted. Therefore, we further tested whether the 
expected Euclidean length of synaptic connections that were most 
likely to be deleted correlated with the current average activity in 
the network. However, an activity-dependent effect on synapse 
deletion was not observed (Figures 9C,D). 

3.3. HOMEOSTASIS DECREASES THE BETWEENNESS CENTRALITY IN 
SMALL-WORLD NETWORKS 

Betweenness centrality (or inbetweenness) is a measure for the 
importance of nodes in a network. Compared with random net- 
works, small-world networks without homeostatic synapse for- 
mation had a relatively high betweenness centrality (Figure 10A). 
Small-world networks with homeostasis, by contrast, revealed 
a pronounced decrease in betweenness centrality over develop- 
mental time. Values peaked in the beginning of development 
when neurons first connected to their nearest neighbors only. 
However, right after the peak, homeostatic synapse formation 
generated topologies with values for betweenness centrality that 
were lower than in non-homeostatic networks. Betweenness 
centrality is high in non-homeostatic small-world networks 
because synapse formation is only determined by the kernel func- 
tion (Equation 5), which often caused that the same cell pairs 
were connected repeatedly. The consequence is a limited num- 
ber of shortest paths between any pairs of nodes in the network. 
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FIGURE 9 | The spatial distribution of newly formed synapses change in 
dependence on the calcium concentration. In (A,B), each dot represents the 
Euclidean distances between those neurons that are most likely to form a 
synaptic connection with each other at this update in connectivity. For this, we 
took at every update in connectivity in which vacant synaptic elements were 
available the Euclidean distance of the connection from neuron j to / for which 
pform (Equation 7) W as maximal. In (CD), each dot represents the length of that 
connection (again in terms of the Euclidean distance between the connected 
neurons) for which synapse deletion was most likely, i.e., Pfj 1 (Equation 6) was 
maximal for every update in connectivity in which synapses had to be deleted. 
In (A,C), we plotted the Euclidean distances for synapse formation and deletion 



By contrast, in homeostatic networks, synapse formation depends 
on the availability of synaptic elements, which can force synapses 
to be formed that are less likely according to the kernel. This 
may increase the number of alternative paths between two neu- 
rons and therefore also increase the number of multiple shortest 
paths through the network. Hence, betweenness centrality quickly 
decreased over time. 



over time. The black curve (right y-axis) indicates the course of the calcium 
concentration [Ca 2 +]. In (B,D), we plotted synapse formation and deletion in 
dependence on [Ca 2+ ] . The color code in all panels indicates the density of the 
dots in the diagrams, with blue and red representing low and high densities of 
dots, respectively. The figure essentially shows that before calcium reaches the 
homeostatic set-point e, the distribution for synapse formation is rather 
Gaussian, following the Kernel function K (Equation 5). The distribution 
becomes random and scattered, with increased Euclidean distances, when 
calcium is at the set-point. The stripes in the distribution arise from the fact that 
not all Euclidean distances are possible due to the grid layout of the network. 
There is no change in the distribution for synapse deletion. 



Because the variety of shortest paths in non-homeostatic net- 
works is limited, betweenness centrality reached a stable plateau 
in these networks. Interestingly, the betweenness centrality in 
homeostatic networks initially also converged toward a quasi- 
stable level. However, as soon as activities approached the home- 
ostatic set-point e, betweenness centrality strongly decreased. 
Over time, the rate of decrease slowed down. The decrease in 
betweenness centrality precisely coincided with the increase in 
Euclidean lengths of synaptic connections. Because the kernel 
function favors short- over long-range connections and therefore 
initially creates networks with high betweenness centrality, we 
may conclude that the formation of longer-range connections (in 
an Euclidean sense) lead to a decreasing betweenness centrality. 
Any new long-range connection not present in the network before 
creates new shortest paths, which in turn decreases betweenness 
centrality. By contrast, homeostasis did not affect the course of 
betweenness centrality in random networks (Figure 10B). 

3.4. SMALL-WORLD NETWORKS BECOME MORE EFFICIENT BY 
HOMEOSTASIS 

Efficient information transmission is probably the most-desired 
property in computational networks. Small-world networks are 
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FIGURE 10 | Changing betweenness centrality over time. (A) Small-world 
networks with (orange) and without (gray) homeostasis in electrical activity. 
(B) Random networks with (blue) and without (gray) homeostasis in electrical 
activity. Means over five simulations per scenario. Shadings of the curves 
indicate standard deviations. 
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very efficient because they combine a high clustering coefficient 
with a short characteristic path length. Nevertheless, their effi- 
ciency is still markedly lower than that of random networks. 
Homeostatic small-world networks, however, generated efficiency 
levels that during the whole course of development exceeded the 
levels in non-homeostatic small-world networks (Figure 11 A). 
Remarkably, at the time when average activities reached the 
homeostatic set-point e, efficiency levels of homeostatic net- 
works further increased and almost reached the levels in random 
networks (Figure 11B). Consequently, favoring more distant con- 
nections in combination with homeostasis in electrical activity 
led to a more efficient network topology than achieved without 
homeostasis. 

4. DISCUSSION 

We have shown that network formation favoring short-range over 
long-range connections produced networks with a pronounced 
small-world structure. Networks with homeostasis in electrical 
activity developed a weaker small-world structure in favor of 
more efficient wiring of connections. Global efficiency particu- 
larly increased when network activity reached the homeostatic 
set-point. Increased global efficiency was caused by the fact that 
homeostasis favored longer-ranged connections, which affected 
clustering as well as characteristic path length. Thus, network 
topology continued to change even after the network had reached 
a homeostatic equilibrium in electrical activity. 

Adding more long-range connections to a small-world net- 
work makes the network more efficient but also more ran- 
dom. This is apparent in our simulations, too, by a decreas- 
ing clustering coefficient and a decreasing characteristic path 
length. Nevertheless, the small-world property of the net- 
works was preserved throughout the whole course of devel- 
opment, and the decrease in clustering coefficient was slow- 
ing down over time. However, networks would most likely 
turn into random networks if rewiring continued indefinitely. 
Consequently, there seems to be a trade-off in network devel- 
opment between high clustering and strong small-worldness on 
the one hand and more randomness and higher efficiency on 
the other hand. The latter particularly arises when networks con- 
tinue to rewire their circuitry when they are in a homeostatic 
equilibrium. 



With additional long-range connections betweenness central- 
ity decreases. Networks with low betweenness centrality are more 
robust against lesions because all nodes are equally important. 
By contrast, networks with high betweenness centrality are very 
vulnerable to lesions. If neurons that are part of many short- 
est paths are lost, the characteristic path length will immediately 
increase, with a significant impact on information transmission. 
It is remarkable that a self-organizing process that forms net- 
works by striving toward homeostasis in electrical activity as a side 
effect produces topologies with lower betweenness centrality that 
contribute to a higher robustness against lesions. 

Key to the homeostasis-driven change in topology is the 
increasing Euclidean length of connections. Since we did not 
observe an increase in connection length in networks with- 
out homeostasis, we concluded that the increase in connection 
lengths was caused by homeostasis in electrical activity. To create 
networks without homeostasis, we took the number of synapses 
from the homeostatic network at every update in connectiv- 
ity and distributed them randomly under the same kernel in 
the non-homeostatic network. There are, in fact, other ways to 
create networks without homeostasis. We additionally built non- 
homeostatic networks by initially giving all model neurons a fixed 
number of vacant synaptic elements and then running updates 
in connectivity until no more synapses could be formed. Also in 
this scenario we did not see an increase in connection length. In 
a third scenario, we added a few vacant synaptic elements to all 
neurons before every update in connectivity. No matter how few 
vacant elements we added and how long we ran the network, we 
did not obtain more long-range connections than expected from 
the kernel. We only observed more long-range connections when 
we slowly added vacant elements at a few randomly selected neu- 
rons at a time (i.e., spatial sparseness in the formation of synaptic 
elements). From these observations we concluded that the contri- 
bution of homeostasis is not only to limit the number of available 
vacant synaptic elements but also to generate a certain sparseness 
in the formation of vacant elements. In fact, as long as the neurons 
were far away from the homeostatic set-point, synaptic element 
formation was not sparse at all as all neurons added elements 
roughly at the same time at equal rates. Only when the neurons 
approached the homeostatic set-point did sparseness in synap- 
tic element formation arise. Homeostasis creates this sparseness 
because balancing the activity of one neuron immediately affects 
the activity in other neurons, which in turn may be driven fur- 
ther away from the homeostatic set-point and then start forming 
new vacant elements. The presence of inhibitory neurons would 
further reinforce this process. 

Homeostasis in electrical activity is one way by which a local 
process can give rise to a change in global topology. Another com- 
parable mechanisms was provided by Kaiser et al. (2009). In their 
model study, they showed that a simple axonal growth process 
can generate the experimentally observed exponential decrease 
in number of connections with increasing connection length. As 
a result of the growth process, most connections become short- 
range, although long-range connections also arise, but in lower 
numbers. The idea of their model is that axons grow out until 
they hit a postsynaptic target. The capacity of model neurons to 
receive connections is limited, and if nearby target neurons are 
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FIGURE 11 | Global efficiency changes over time due to network 
development. (A) Small-world networks with (orange) and without (gray) 
homeostasis in electrical activity. (B) Random networks with (blue) and 
without (gray) homeostasis in electrical activity. Means over five simulations 
per scenario. Shadings of the curves indicate standard deviations. 



Frontiers in Synaptic Neuroscience 



www.frontiersin.org 



April 2014 | Volume 6 | Article 7 | 11 



Butz et al. 



Homeostatic small-world networks 



completely occupied by incoming connections, axons continue to 
grow out until they hit a vacant target. Hence, Euclidean connec- 
tion length increases over developmental time. In this model, the 
spatial growth process in combination with a hard boundary on 
the number of incoming synapses per model neuron generates the 
increase in connection lengths. 

There are striking similarities between the topology of our 
model networks and that of developing dissociated cell cultures. 
It is well known that cultured neuronal networks can form small- 
world topology (Bettencourt et al, 2007; Yu et al, 2008; Gerhard 
et al, 2011). Downes et al. (2012) reported an increase in clus- 
tering coefficient of dissociated cell cultures between 14 days in 
vitro (DIV) and 28 DIV and a subsequent drop until 35 DIV, a 
course of development that is comparable to the course of devel- 
opment in our model networks with homeostasis. Between 14 
and 35 DIV, the mean shortest path length did not change signif- 
icantly, only showing a slight drop around 28 DIV. Consequently, 
small- worldness reached its maximum around 28 DIV. Moreover, 
the experimental data indicated the presence of longer synaptic 
connections from 28 DIV onwards that had not been not present 
at 2 1 DIV. From our previous studies we know that dissociated 
cell cultures reach homeostasis around this time (Tetzlaff et al., 
2010). Therefore, we may hypothesize that the increase in con- 
nection length in dissociated cell cultures may be due to neurons 
reaching a homeostatic equilibrium in electrical activity. Other 
synaptic plasticity mechanisms not currently incorporated in our 
model may of course also have contributed to network formation 
in developing cell cultures. 

On a macroscopic scale, functional imaging data reveal a devel- 
opment of small-world topology that also has interesting similar- 
ities with the self-organizing network formation in our model. 
The infant human brain has small-world properties already at the 
third post-natal week (Fransson et al., 2011). During the follow- 
ing 2 years, network topology undergoes a significant refinement: 
brain networks increase their small-worldness, global efficiency 
and number of long-distance connections (Gao et al., 2011). 
Although our network model shows only an initial transient 
increase in small-worldness, it may offer a simple explanation for 
the sudden increase in number of long-range connections and 
the associated increase in global efficiency. Could it be that even 
in the human brain, neurons establishing a homeostatic equilib- 
rium in electrical activity produce — as an emergent property of 
the homeostatic growth process — more long-range connections? 
Remarkably, the increase in number of long-range connections 
occurs not during a genetically-encoded formation of an initial 
embryonic layout of projections but during the post-natal critical 
period (Gao et al., 2011), during which neurons are highly sensi- 
tive to afferent input. In general, the importance of critical periods 
is to balance excitatory and inhibitory circuits and to establish 
homeostasis in neuronal electrical activity (Hensch, 2005; Butz 
et al, 2009b). Considering that long-range connections arise in 
local as well as global networks, our study raises the interesting 
hypothesis that homeostasis in electrical activity may be the driv- 
ing force for the formation of long-range connections on both a 
microscopic and a macroscopic scale. 

Homeostasis in electrical activity is a ubiquitous principle in 
the nervous system (Wolff and Wagner, 1983; Ramakers et al., 



1991; Abbott and Nelson, 2000) and a variety of plasticity mech- 
anisms can act homeostatically. Scaling of synaptic strengths 
(Turrigiano, 1999), for example, has been reported as a mech- 
anism acting at existing synapses to stabilize postsynaptic firing 
in cortical, hippocampal and spinal cord networks (Lissin et al., 
1998; O'Brien et al., 1998; Turrigiano and Nelson, 1998). Even 
in the mature brain, not only the strength of synapses but also 
the formation of new synapses can contribute to the stabilization 
of neuronal activity, for example after focal retinal lesions (Butz 
and van Ooyen, 2013). In developing dissociated cell cultures, we 
showed that homeostasis in electrical activity may be a precon- 
dition for the emergence of self-organized criticality in neuronal 
firing (Tetzlaff et al., 2010). In another study, we showed that 
homeostasis in electrical activity can regulate the synaptic embed- 
ding of newly formed neurons in the mature hippocampal dentate 
gyrus (adult neurogenesis) and can account for the experimen- 
tally observed counter-intuitive inverse relationship between cell 
proliferation rate and synaptogenesis (Butz et al., 2008). 

In summary, we conclude that homeostatic regulation of elec- 
trical activity together with simple distance-dependent forma- 
tion of connections is capable of creating, in a self-organizing 
manner, neuronal networks that are more robust and more 
efficient than networks grown without homeostatic regulation. 
Strikingly, the growth process revealed features of developing 
topologies that are also observed in dissociated cell cultures and 
infant human brains. The formation of network topology by a 
self-organizing, local growth process may also be relevant for 
automatically generating the connectivity structure of large-scale 
neuronal networks (Potjans and Diesmann, 2014) that are cur- 
rently studied in enterprizes such as the Human Brain Project 
(www.humanbrainproject.eu). 
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